#   LibAut
#       Figure 3 
#   Juraj Medzihorsky
#   2021-08-10


options(StringsAsFactors=F)
gr <- function(i) ((1+sqrt(5))/2)^i

#   sessionInfo()
##  R version 3.6.3 (2020-02-29)
##  Platform: x86_64-pc-linux-gnu (64-bit)
##  Running under: Ubuntu 16.04.7 LTS
##
##  Matrix products: default
##  BLAS/LAPACK: /opt/OpenBLAS/lib/libopenblas_haswellp-r0.3.9.so
##
##  locale:
##   [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
##   [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_GB.UTF-8       LC_NAME=C
##   [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
##  attached base packages:
##  [1] stats     graphics  grDevices utils     datasets  methods   base
##
##  other attached packages:
##  [1] nvimcom_0.9-28 colorout_1.2-0
##
##  loaded via a namespace (and not attached):
##  [1] compiler_3.6.3 tools_3.6.3


#-------------------------------------------------------------------------------
#   load the data

code_dir <- getwd()
data_dir <- gsub('code$', 'data', code_dir)
figs_dir <- gsub('code$', 'figures', code_dir)


setwd(data_dir)
e <- read.csv('eps.csv')    
setwd(figs_dir)

e$dem_ep_id <- as.character(e$dem_ep_id)

#   mark the libaut episodes
e$cy <- paste0(e$country_text_id, '_', e$year)
e$libaut <- ifelse((e$dem_ep_prch%in%1) & (e$dem_ep_ptr%in%0), 1, 0)
e$libaut_id <- rep('0000', nrow(e))
e$libaut_id[e$libaut%in%1] <- e$dem_ep_id[e$libaut%in%1]

with(e, table(dem_ep_prch, dem_ep_ptr, libaut))


#-------------------------------------------------------------------------------
#   prep and plot

#   dem_ep_prch
#       0   not in an ep with dem transition potential
#       1   in an ep with potential for a dem transition

#   dem_ep_ptr: 
#       0   not in an a dem transition episode
#       1   in an dem transition episode

#   only episodes that start with sub_dem
good_fil <-
    (e$dem_ep_prch%in%1) & (e$dem_ep_ptr%in%0) &
    (!is.na(e$dem_ep_start_year)) &
    (e$year==e$dem_ep_start_year)  
good_eps <- e$dem_ep_id[good_fil]



E <- e
E$new_dem_id <- E$libaut_id

E$chunk <- rep(NA, nrow(E))

#   walk line by line, add chunk indicator
#   this assumes the data is ordered by country 
#   and within country by year
E$chunk[1] <- 1
for (i in 2:nrow(E)) {
    cond1 <- E$year[i]!=(E$year[i-1]+1)                     #   year jump
    cond2 <- E$country_text_id[i]!=E$country_text_id[i-1]   #   different country
    cond3 <- E$new_dem_id[i]!=E$new_dem_id[i-1]             #   ep start/end (incl preep)
    if (cond1|cond2|cond3) {
        E$chunk[i] <- E$chunk[i-1]+1
    } else {
        E$chunk[i] <- E$chunk[i-1]
    } 
}


#   get tables for barplots
#       ctab    continuing
#       otab    originating
atab <- with(E, table(year, dem_ep_outcome))
ctab <- data.frame('year'=rownames(atab),
                   'nope'=atab[,1], 
                   'dem'=rowSums(atab[,c(2,6)]),
                   'aut'=rowSums(atab[,c(3:5)]), 
                   'cen'=atab[,7])


btab <- with(e[(!is.na(e$dem_ep_start_year))&(e$year==e$dem_ep_start_year),], 
             table(dem_ep_start_year, dem_ep_outcome))
otab <- data.frame('year'=rownames(btab),
                   'dem'=rowSums(btab[,c(1,5)]),
                   'aut'=rowSums(btab[,c(2:4)]), 
                   'cen'=btab[,6])


ec <- split(E, f=E$chunk)


#           blue        orange    green
eks <- c( '#56B4E9', '#E69F00', '#009E73')

#   #   DEM_EP_OUTCOME
#   1   0– There is no democratization episode during this year
#   2   1– Democratic transition 
#   3   2– Preempted democratic transition 
#   4   3– Stabilized electoral autocracy 
#   5   4– Reverted liberalization 
#   6   5– Democratic deepening 
#   7   6– Censored

#   #   recode outcomes
#   #   epy: suc 1 fail 0
#   epframe$epy <- as.numeric(rep(NA,nrow(epframe)))
#   #   success: 1, 5
#   epframe$epy[epframe$final_outcome%in%c(1,5)] <- 1
#   #   fail: 2,3,4
#   epframe$epy[epframe$final_outcome%in%c(2,3,4)] <- 0
#   #   censored: 6
#   epframe$epy[epframe$final_outcome%in%6] <- NA


#   function to plot a single chunk
plotchunk <-
    function(x, a=0.5) 
    {
        nr <- nrow(x)
        if (x$new_dem_id[1]%in%good_eps) {
            if (x$dem_ep_outcome[nr]%in%0) {
                tkl <- grey(0.8, 0.3)
                twd <- gr(-1)
            } else if (x$dem_ep_outcome[nr]%in%c(1,5)) {
                tkl <- paste0('#56B4E9', as.hexmode(floor(a*255)))
                twd <- gr(0) 
            } else if (x$dem_ep_outcome[nr]%in%c(2,3,4)) {
                tkl <- paste0('#E69F00', as.hexmode(floor(a*255)))
                twd <- gr(0)
            } else if (x$dem_ep_outcome[nr]%in%6) {
                tkl <- paste0('#009E73', as.hexmode(floor(a*255)))
                twd <- gr(0)
            } 
        } else {
            tkl <- grey(0.8, 0.3)
            twd <- gr(-1)
        }
        #   extend lines after
        tx <- x$year
        ty <- x$v2x_polyarchy
        #
        #   if a non-episode starts and an episode ended, bridge the gap
        co1 <- x$new_dem_id[1]%in%'0000'
        co2a <- E$country_text_id%in%x$country_text_id[1]
        co2b <- E$year%in%(x$year[1]-1)
        co2 <- sum(co2a&co2b)==1
        if (co2) {
            fx <- x$year[1]-1
            fy <- E$v2x_polyarchy[co2a&co2b]
            tx <- c(fx, tx)
            ty <- c(fy, ty)
        }
        #   same, but for the end
        co2c <- E$year%in%(x$year[nrow(x)]+1)
        co3 <- sum(co2a&co2c)==1
        if (co1&co3) {
            lx <- x$year[nrow(x)]+1
            ly <- E$v2x_polyarchy[co2a&co2c]
            tx <- c(tx, lx)
            ty <- c(ty, ly)
        }
        #
        lines(tx, ty, col=tkl, lwd=twd) 
        #   start end dots
        if (x$dem_ep_outcome[nr]%in%c(1:5)) {
            l <- nrow(x)
            l <- length(tx)
            points(tx[c(1,l)], ty[c(1,l)], pch=16, cex=0.38, col=tkl)
        } else if (x$dem_ep_outcome[1]%in%c(6)) {
            points(tx[1], ty[1], pch=16, cex=0.38, col=tkl)
        }
    }


#   now plot

xl <- c(1900, 2020) + c(-1,1)*c(1.618, 0.618)
yl <- c(0.38,1)
yl_true <- c(-1-1.618^-7,1)
y0 <- seq(0,1,by=0.5)
y1 <- seq(0, 10, by=10)
y2 <- seq(0, 200, by=100)
png('Figure_3_raw.png', 1300*3, 900*3, res=300)
par(mfrow=c(1,1), mar=c(3,1,1,3), xpd=F, xaxs='i', yaxs='i', xpd=T)
plot(0,0,type='n',frame=F,xaxt='n',yaxt='n',xlab='',ylab='',main='',
     xlim=xl,ylim=yl_true)
    x0 <- seq(1900,2020,by=20)
    axis(1,x0,rep('',length(x0)),lwd=0.38, tck=-0.382*1e-2)
    axis(1,x0,x0,lwd=0,line=-0.618)     
    text(1900-0.5, 1-1.618^-6, 'Liberalization Episodes and Electoral Democracy', pos=4, cex=1.235, offset=0)
    text(1900-0.5, 0-1.618^-6, 'Episode Onsets', pos=4, cex=1.235, offset=0)
    text(1900-0.5, -0.38-1.618^-6, 'Countries in Episodes', pos=4, cex=1.235, offset=0)
    axis(4,y0,rep('',3),lwd=0.38,las=2, tck=-0.382*1e-2)
    axis(4,y0,y0,lwd=0.0,las=2, line=-0.618)
invisible(sapply(ec, plotchunk))
#
    for (i in 1:nrow(otab)) {
        anchor <- -0.19
        stretch <- abs(anchor)/18.525
        tx <- as.numeric(rownames(otab)[i])+0.0+c(-1,1)*1.618^-2
        ty3 <- anchor + c(0, c(otab[i,2]+otab[i,4])*stretch)
        ty1 <- anchor + c(0, otab[i,2]*stretch)
        polygon(c(tx,rev(tx)), c(ty3[1],ty3[1],ty3[2],ty3[2]), lwd=gr(-1), col=eks[3], border='white')
        polygon(c(tx,rev(tx)), c(ty1[1],ty1[1],ty1[2],ty1[2]), lwd=gr(-1), col=eks[1], border='white')
    }
    for (i in 1:nrow(otab)) {
        anchor <- -0.19
        stretch <- abs(anchor)/18.525
        tx <- as.numeric(rownames(otab)[i])+0.0+c(-1,1)*1.618^-2
        ty <- anchor + c(-otab[i,3]*stretch, 0)
        polygon(c(tx,rev(tx)), c(ty[1],ty[1],ty[2],ty[2]), lwd=gr(-1), col=eks[2], border='white')
    }
    abline(h=-0.38/2, col='white')
    axis(4,anchor +  y1*stretch, rep('',2), las=2, lwd=0.38, tck=-0.382*1e-2) 
    axis(4,anchor + -y1*stretch, rep('',2), las=2, lwd=0.38, tck=-0.382*1e-2) 
    axis(4,anchor +  y1*stretch, y1, las=2, lwd=0.0, line=-0.618) 
    axis(4,anchor + -y1*stretch, y1, las=2, lwd=0.0, line=-0.618) 
#
for (i in 1:nrow(ctab)) {
    anchor <- -1
    stretch <- 0.62/200
    tx <- as.numeric(rownames(ctab)[i])+0.0+c(-1,1)*1.618^-2
    ty0 <- anchor + c(0, rowSums(ctab[,-1])[i]*stretch)
    ty2 <- anchor + c(0, sum(ctab[i,3:5])*stretch)
    ty3 <- anchor + c(0, sum(ctab[i,c(3,5)])*stretch)
    ty1 <- anchor + c(0, sum(ctab[i,3])*stretch)
    polygon(c(tx,rev(tx)), c(ty0[1],ty0[1],ty0[2],ty0[2]), lwd=gr(-1), col=grey(0.8, 0.5), border='white')
    polygon(c(tx,rev(tx)), c(ty2[1],ty2[1],ty2[2],ty2[2]), lwd=gr(-1), col=eks[2], border='white')
    polygon(c(tx,rev(tx)), c(ty3[1],ty3[1],ty3[2],ty3[2]), lwd=gr(-1), col=eks[3], border='white')
    polygon(c(tx,rev(tx)), c(ty1[1],ty1[1],ty1[2],ty1[2]), lwd=gr(-1), col=eks[1], border='white')
}
axis(4,anchor + y2*stretch, rep('',3), las=2, lwd=0.38, tck=-0.382*1e-2) 
axis(4,anchor + y2*stretch, y2, las=2, lwd=0.0, line=-0.618) 
#
par(xpd=T)
legend('topleft', inset=c(0,+0.0553), horiz=F,
       pch=rep(15,3), col=c(eks[c(1,3,2)], grey(0.8,0.5)),
       pt.cex=rep(1.618^1,4), bty='n', cex=rep(1-0.235,4),
       legend=c('Succesful', 'Censored', 'Failed', 'None'))
par(xpd=F)
#
dev.off()

#   cropping the figure 
#   -   doesn't affect figure content
#   -   requires additional software installed
system('convert -trim Figure_3_raw.png Figure_3.png')

#   SCRIPT END